Transport properties for driven granular fluids in situations close to steady 

homogeneous states 
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The transport coefficients of a granular fluid driven by a stochastic bath with friction are obtained 
by solving the inelastic Enskog kinetic equation from the Chapman-Enskog method. The heat and 
momentum fluxes as well as the cooling rate are determined to first order in the deviations of 
the hydrodynamic field gradients from their values in the homogeneous steady state. Since the 
collisional cooling cannot be compensated locally for the heat produced by the external driving 
force, the reference distribution (zeroth-order approximation) depends on time through its 
dependence on temperature. This fact gives rise to conceptual and practical difficulties not present 
in the undriven case. On the other hand, to simplify the analysis and given that we are interested 
in computing transport in the first order of deviations from the reference state, the steady-state 
conditions are considered to get explicit forms for the transport coefficients and the cooling rate. A 
comparison with recent molecular dynamics simulations for driven granular fluids shows an excellent 
agreement for the kinematic viscosity although some discrepancies are observed for the longitudinal 
viscosity and the thermal diffusivity at large densities. Finally, a linear stability analysis of the 
hydrodynamic equations with respect to the homogeneous steady state is performed. As expected, 
no instabilities are found thanks to the presence of the external bath. 
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nal modes in a driven granular fluid. Given that the 
expressions for the transport coefficients were not known 
in this driven problem, the simulation data were com- 
pared with their corresponding elastic forms. Thus, it 
would be desirable to provide simulators with the appro- 
priate theoretical tools to work when studying problems 
in granular fluids driven by thermostats. 

The aim of this paper is to determine the transport co- 
efficients of a dense driven granular gas of inelastic hard 
spheres in the framework of the Enskog kinetic equation. 
As in the undriven case the transport coefficients 
are obtained by solving the Enskog equation by means 
of the Chapman-Enskog expansion 14j around a certain 
reference state /^°^ (zeroth-order approximation). While 
in the undriven case the distribution / W is chosen to be 
the local version of the homogeneous cooling state (HCS), 
there is some flexibility in the choice of /(°) for a driven 
gas. For simplicity, one possibility is to take a local ther- 
mostat such that the distribution /(°) is still stationary 
at any point of the system. This was the choice assumed 
in previous works (la , Il6j to compute the transport coef- 
ficients of a heated granular gas. On the other hand, for 
general small deviations from the steady reference state, 
the zeroth-order distribution /(°) is not in general a sta- 
tionary distribution since the collisional cooling cannot 
be compensated locally by the heat injected by the driv- 
ing force. This fact introduces additional difficulties not 
present in previous studies [15|, Il6| . In this paper, we 
will adopt this point of view and will consider this kind 
of thermostat that seems to be closer to the one used in 
computer simulations. 

The underlying motivation of such a study is twofold. 
* I vicenteg@unex.es|ptp://www.eweb. unex.es/eweb/fisteor/vicente/ 1 On of the one hand, given that most of the computer sim- 
t |fvega@unex.e s| |http://www.eweb.unex.es/eweb/fisteor/fran/ 1 ulations [][ are performed by driving the fluid by means 



I. INTRODUCTION 

It is well established that when a granular material 
is externally excited (rapid flow conditions) the mo- 
tion of grains resembles the random motion of atoms or 
molecules in an ordinary gas. In these conditions, kinetic 
theory together with numerical simulations are the best 
tools to describe the behavior of granular flows. How- 
ever, in contrast to ordinary fluids, the collisions between 
grains are inelastic and so, one has to feed energy into the 
system to achieve a steady non-equilibrium state. This 
can done either by driving through the boundaries, for 
example, shearing the system or vibrating its walls [1|, 
or alternatively by bulk driving, as in air-fluidized beds 
[U, El- O n the other hand, this way of supplying en- 
ergy causes in most of the cases strong spatial gradients 
in the system. To avoid the difficulties associated with 
inhomogeneous states, it is quite usual in computer sim- 
ulations to homogeneously heat the system by the action 
of an external driving force [^-l7|. Borrowing a terminol- 
ogy often used in nonequilibrium molecular dynamics of 
ordinary fluids [8j, this type of external forces are usu- 
ally called "thermostats". Nevertheless, in spite of its 
practical importance, the understanding of the effect of 
the external driving force on the dynamical properties 
of the system (such as the transport coefficients) is still 
not completely understood (9l4ll|. In particular, recent 
molecular dynamics (MD) simulations [y, [7fl have com- 
puted some transport coefficients by measuring the static 
and dynamical structure factors for shear and longitudi- 
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of a thermostat, it is important to know the effect of 
thermostat on the transport properties. In this sense, 
the question arises then as to whether, and if so to what 
extent, the conclusions drawn in previous works [l5l Il6j 
may be altered when the time-dependent distribution 
/(°)(r, v,i) is considered in the Chapman-Enskog expan- 
sion instead of the steady distribution. On the other 
hand, it is also of interest to evaluate the transport co- 
efficients in a situation where a direct comparison with 
computer simulations is more likely to occur. Here, in 
the same way as the simulations carried out in Ref. Q, 
we have assumed that the fluid is driven by a stochastic 
bath with friction. This allows us to compare directly 
our theoretical predictions for the kinematic and longi- 
tudinal viscosities and the thermal diffusivity with those 
obtained from MD simulations Q- 

The plan of the paper is as follows. In Sec. [Til the 
Enskog kinetic equation for a granular gas fluidized by 
an external force is introduced. Then, before considering 
inhomogeneous situations, Sec. IIIII analyzes the steady 
homogeneous state. As in the HCS [IT}, a scaling solu- 
tion ip s is proposed at the steady state. However, the new 
feature is that the dependence of the reduced distribution 
ip s on the granular temperature occurs through two di- 
mensionless parameters (dimensionless velocity and the 
reduced noise strength) instead of a single parameter (di- 
mensionless velocity) as in the HCS. Once the steady ho- 
mogeneous distribution is well characterized, Sec. [IV] ad- 
dresses the Chapman-Enskog expansion around the un- 
steady reference distribution (r, v; t). The details of 
the calculations to first order in the spatial gradients are 
displayed along several Appendices and the explicit ex- 
pressions of the transport coefficients and the cooling rate 
are displayed in Sec. |Vl In reduced form, they are given 
in terms of the volume fraction, the coefficient of restitu- 
tion and the parameters of the thermostat. A comparison 
with recent MD simulations Q for hard disks is done in 
Sec. I VII showing an excellent agreement in the case of the 
kinematic viscosity and some discrepancies for the longi- 
tudinal viscosity and the thermal diffusivity, specially at 
large solid fractions. Section IVIII is devoted to the linear 
stability analysis around the steady homogeneous state. 
As in previous studies 0, 0] based on the elastic forms 
of the Enskog transport coefficients, our results also indi- 
cate that the steady homogeneous state is linearly stable. 
The paper is closed in Sec. IVIIll with a discussion of the 
results derived here. 



II. ENSKOG KINETIC THEORY FOR DRIVEN 
SYSTEMS 

We consider a system of inelastic hard spheres in d di- 
mensions with mass to and diameter a. Collisions are 
characterized by a (constant) coefficient of normal resti- 
tution < a < 1, with a = 1 in the elastic limit. As 
said in the Introduction, in order to maintain a station- 
ary nuidized state, the granular fluid is driven by means 



of an external force or thermostat that acts locally on 
each particle [HI]. This is a quite usual choice in com- 
puter simulations In these conditions, the equation 
of motion for a particle i with velocity can be written 
as 



mv, = Ff (i) 



TiCOll 



(1) 



where F* h is the thermostat force and F™ 11 is the force 
due to inelastic collisions. Here, we assume that F' h is 
composed by two different terms: (i) a stochastic force 
where the particles are randomly kicked between colli- 
sions [l8| and (ii) a viscous drag force which mimics 
the interaction of the particles with and effective viscous 
"bath" at temperature Tb. More explicitly, F* h is given 



by 



Pj h (t) = -7bv < (t)+Ff(*), 



(2) 



where 7b is a drag coefficient that defines the characteris- 
tic interaction time with the external bath, t^ 1 = j^/m. 
As usual, the stochastic force F f* is assumed to have the 
form of a Gaussian white noise [18| : 

(Ff (t)) = 0, (Ff (t)Ff (0) = rn^lkAt ~ O*, (3) 

where 1 is the d x d unit matrix and £j~ represents the 
strength of the correlation. The forcing term in the 
Enskog equation associated to Ff* is represented by a 
Fokker-Planck operator [ijj of the form — ^ 2 d 2 /dv 2 . 
One of the advantages of using the model ^ instead of 
other kind of thermostats is that the temperature of the 
thermostat Tb (different from the kinetic temperature of 
the fluid T < Tb) is always well defined. In particular, 
this kind of thermostat is able to equilibrate the system 
when collisions are elastic. Moreover, a similar exter- 
nal driving force to that of Eq. @ has been recently 
proposed to model the effect of the interstitial fluid on 
grains in monodisperse gas-solid suspensions 19]. 

Thus, the corresponding Enskog kinetic equation for 
the one-particle velocity distribution function /(r,v, t) 
reads 

7b ,, d 7b d 

dtf + v-Vf U • — / — • V/ 

to oV m oV 

-^b^/ = ^[r,v|/,/], (4) 

V = v — U is the peculiar velocity, U is the mean flow 
velocity, and 

j E [r, Vl |/,f] = ° d ~ 1 J rfv 2 J dae(a- Sl2 )(a- gl2 ) 

x [a- 2 X {r, r - cr)/(r, v[ ; t)f(r - <x, v 2 ; <) 
- X (r,r + < r)/(r,v 1 ;t)/(r + < T,v 2 ;t)] (5) 

is the Enskog collision operator. Here, d is the dimen- 
sionality of the system (d = 2 for disks and d = 3 for 
spheres), er = <7<x, <x being a unit vector, O is the Heavi- 
side step function, and gi2 = vi — v 2 . The primes on the 
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velocities in Eq. (J5)) denote the initial values {v^, v 2 } that 
lead to {v l7 v 2 } following a binary collision, v' x = v x — 
| (1 + a" 1 ) (ct ■ gi 2 )S, v 2 = v 2 + | (1 + a" 1 ) (ct ■ gl2 )ff. 
In addition, x[r, r + er|{n(t)] is the equilibrium pair corre- 
lation function at contact as a functional of the nonequi- 
librium density field n(r, t) defined by 



n(r,t) = J efv/(r,v,t). 



(0) 



The macroscopic balance equations for the system are 
obtained when one multiplies the Enskog equation <(4]) 
by {1, mv, mv 2 } and integrates over velocity. After some 
algebra one gets [HI, [ll| 



D t n + nV ■ U = 



D t V = -/3 -1 V • P 



7b 



U 



(7) 



(8) 



D t T+ — (V-q+P: VU) 

an 



2T 



Tb + Kb-CT. (9) 



In the above equations, D t = dt + U • V is the material 
derivative and p — mn is the mass density. The cooling 
rate £ is proportional to 1 — a 2 and is due to dissipativc 
collisions. The pressure tensor P(r,t) and the heat flux 
q(r, t) have both kinetic and collisional transfer contri- 
butions, i.e., P = P fe + P c and q = q fc + q c . The kinetic 
contributions are given by 

P fc = J dvmVV/(r,v,i), q fe = J dvyy 2 V/(r,v,i), 

_ (io) 

and the collisional transfer contributions are |12| 
P c = i±^ma d J d Vl Jdv 2 J rfSe(5-g 12 )(?-g 12 ) 2 



; / dx / (2) [r - xa,r + (1 - x)cr, Vi, v 2 ; t] 
Jo 



(11) 



1 + a 



ma d J dvi J dv 2 J da 9(<T • gi 2 )(<r • g i2 ) 2 

x(Gi 2 -<r)er / [r — xcr, r + (1 — x)(t, vi, v 2 ; t] . 

Jo 

(12) 

Here, Gi 2 = |(Vi + V 2 ) is the velocity of center of mass 
and 

/ (2) (ri,r 2 ,vi,v 2 ,f) = x(ri,r 2 |n(t))/(ri,vi,t)/(r 2 ,v 2 ,t). 

(13) 



Finally, the cooling rate is given by 

-mcr^ 1 f dvi / dv 2 I da@(a ■ gi 2 ) 



^ 4dnT 

x(5 ; - gl2 ) 3 / (2) (r,r + ^,v 1 ,v 2 ;t). 



(14) 



III. STEADY HOMOGENEOUS STATES 

Before considering inhomogeneous situations, it is 
quite instructive to analyze first the homogeneous state. 
In this situation, the density n(r, t) = n s is constant, the 
granular temperature T(r, t) = T(t) is spatially uniform 
and the mean flow vanishes (U = 0). As a consequence, 
the one-particle distribution function f(v,t) verifies the 
kinetic equation 



dtf 



7b 9 
m dv 



v/-^^/=Jb[/,/], (15) 



where 



Je [/, /] - x^ 1 J rfv 2 J d*e(B ■ gl2 )(5 • gl2 ) 

x [«- 2 /K)/(«a)-/(«i)/(«a)]. (16) 

Here, \ is t ne P a i r correlation function evaluated at the 
(homogeneous) density n s . The collision operator (fi"6j) 
can be recognized as the Boltzmann operator for inelastic 
collisions multiplied by the factor \. The energy balance 
equation reads simply 



dtT 



2T 



-7b + m£, h - C T. 



(17) 



In the hydrodynamic regime, / qualifies as a normal so- 
lution and so, its time dependence only occurs through 
the granular temperature T: 



df 

dtf = TT-ftT 



m 



7b - ~tl + C) T 



9/ 
dT' 



(18) 



Substitution of Eq. (JTSJ) into Eq. JTSJ yields 

— 7b - 777 Sb + C — • V/ 

m J J oT m ov 

1 <9 2 

~ 2^a^ / = Je[/ ' /] ' 

After a transient regime, the gas will achieve a steady 
state characterized by a constant temperature T s . Ac- 
cording to Eq. (|17p , the steady granular temperature T s 
is given by the equation 



(19) 



27b, 
m 



( S T S + —T s = m£b, 



(20) 



where the (steady) cooling rate £ s is defined by Eq. 
(1141) by using the stationary distribution function / s (v). 
Equation ([20)) establishes a relation between the model 
parameters 7b and £ 2 so that only one of the above pa- 
rameters is independent in the steady state. Here, we will 
take the noise strength £^ as the relevant external param- 
eter. By using the relation (|2T))) . in the steady state Eq. 
(IT91 becomes 



v/ s -^b|^/s = Je[/ s ,/ s ]. (21) 



1.9 , m^ 2 d 

— Ct V 7= 

2 s 9v J 2T S dv 
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Equation ([21]) shows clearly that the steady distribu- 
tion f s also depends on the model parameter ^ (apart 
from its dependence on the coefficient of restitution and 
the granular temperature). Thus, although the explicit 
form of f s is not known so far, dimensionless analysis 
requires that / s has the scaled form [2(3] 



f,(v,&)=n a Vo d <p, (c,C) 



(22) 



where (p s is an unknown function of the dimensionless 
parameters 



— . £ 



T s v 



(23) 



Here, vq = \j2T s /m is the thermal speed and I = 
l/(n s cr d_1 ) is the mean free path for hard spheres. Note 
that the dependence of the scaled distribution function 
ip s on the temperature is encoded through two parame- 
ters: the dimensionless velocity c and the (reduced) noise 
strength £* . This scaling differs from the one assumed in 
the case of the HCS [T7J where only the dimensionless 
velocity c is required to characterize the distribution <p a , 
A similar scaling solution to the form (|2"2")l has been re- 
cently proposed [21| for a driven homogeneous granular 
gas before reaching the stationary regime. In terms of 
the (reduced) distribution function <p a , Eq. (|21[) can be 
finally rewritten as 



1~ d 



1,» d 



1 €^><P S = Je[Vs,Vs], (24) 



4" s d 



where we have introduced the dimensionless quantities 



'as — ) 

vo 



T* 



- Je- 



(25) 



Since the cooling rate vanishes for elastic collisions, 
then the solution of Eq. (pM)) is the Maxwellian distribu- 
tion 



<p B (c) 



-d/2 - 

7r ' e 



(26) 



However, if the particles collide inelastically (a < 1), 
C* 7^ and the exact form of (p s (c) is not known. In 
particular, the deviation of <p s (c, £*) from its Maxwellian 
form (|26[) is measured through the kurtosis or fourth- 
cumulant 



a-2.: 



d(d + 2) 



1, 



where 



dc c k ip s (c). 



The steady-state value a2, s of the kurtosis can be deter- 
mined by considering the leading Sonine approximation 
for Q: 



T d/2 



1 + 02,i 



(d + 2)c 2 d(d + 2) 
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FIG. 1. Plot of the reduced granular temperature T s /Tb ver- 
sus the volume fraction (j> for a two-dimensional [d = 2) gran- 
ular fluid and two different values of the coefficient of resti- 
tution: a = 0.8 (solid line) and a = 0.6 (dashed line). The 
symbols are the Monte Carlo simulation results (circles for 
a = 0.8 and triangles for a — 0.6). 




FIG. 2. Plot of the steady fourth cumulant 02, s versus the 
coefficient of restitution a for a two-dimensional (d = 2) gran- 
ular fluid with cj> — 0.25. The line is the theoretical result and 
the symbols are the Monte Carlo simulation results. 



The approximation (|29|) is justified because the coeffi- 
cient a,2,s is expected to be small [17]. With the use of 
the form (|29p and neglecting nonlinear terms in d2. s , the 
dependence of 02 lS on a and £* can be explicitly deter- 
mined. After some algebra, one gets (2(| 



02,1 



16(1 - a)(l - 2a 2 ) 



9 + 24d - a(41 - 8d) + 30(1 - a)a 2 + C d - 



x(i+«) 
(30) 

(27) where C d = 16d(d+2)\/2r(d/2)/7r( d - 1 )/ 2 . In the absence 
of friction (7b = 0), the steady-state condition yields 
C s * = C an d Eq- ([30t agrees with the results obtained 
when the system is only driven by the stochastic ther- 

(28) mostat ljj. Moreover, when_ 



* = 0, we also recover the 
results of the undriven case [17] . 

Once the kurtosis is known, the cooling rate £ s can be 
written in terms of 02, s as 



Cs = - ^ (1 ~ a 2 )x ( 1 + ^a 2 ,s ) n a c7 d - 1 - 



d r. 



m 
( 31 ) 

(29) where the steady granular temperature T s obeys the 
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equation 

m 2 il 2 d " 1 [mx<t>, 3 \ ,„ 

T s = 2k . /_*T(i - a 2 ) 1 + — a 2 B Ti 5 / 2 . 

(32) 

d (33) 



Here. 



T d/2 



" 2 d - 1 dr(f) s 

is the solid volume fraction. Equation (1321) gives the gran- 
ular temperature T s in the non-equilibrium stationary 
state. 

Figure Q] shows the (reduced) steady temperature 
T s /Tb versus the volume fraction <p for two different val- 
ues of the coefficient of restitution a. The theoretical 
results obtained from Eq. (f3"2"|) for hard disks (d = 2) are 
compared with those obtained by numerically solving the 
Enskog-Boltzmann equation from the direct simulation 
Monte Carlo (DSMC) method [H]. As in Ref. 0, the 
fixed parameters of the simulations are m = 1, a = 0.01, 
7b = 1, 7b = 1, and £ 2 = 2. For a two-dimensional 
system, we have chosen the following form for x{4>) [13] : 



X(<t>) 



1-^ 

(1-0) 



(34) 



We observe an excellent agreement between theory and 
simulation in the complete range of values of <p consid- 
ered. As expected, at a given value of the solid fraction, 
the steady granular temperature T s decreases as the gas 
becomes more inelastic. The dependence of the kurto- 
sis ct2. s on a is shown in Fig. [2] for <f> = 0.25, d = 2, 
and the same simulation parameters as in Fig. [TJ It is 
quite apparent that simulation data compare very well 
with the theoretical result (|30|) , even for extreme values 
of dissipation. This good agreement suggests that the 
scaled distribution </5 s (c, £*) can be well represented by 
the leading Sonine approximation (|29[) in the region of 
thermal velocities. In addition, the values of a2, s in the 
driven case are generally smaller than in the undriven 
case El SI. 



adapted to dissipative dynamics. The Chapman-Enskog 
method assumes the existence of a normal solution such 
that all space and time dependence of the distribution 
function occurs through the hydrodynamic fields 



/(r,v,t)=/[v|n(r,t),T(r,t),U(r,t)] 



(35) 



The notation on the right hand side indicates a func- 
tional dependence on the density, temperature and flow 
velocity. For small spatial variations (i.e., low Knudsen 
numbers), this functional dependence can be made local 
in space through an expansion in the gradients of the 
hydrodynamic fields. To generate it, / is written as a 
series expansion in a formal parameter e measuring the 
non-uniformity of the system, 



/ = / (0) +£/ (l) +e 2 /( 2) + 



(36) 



where each factor of e means an implicit gradient of a 
hydrodynamic field. The uniformity parameter e is re- 
lated to the Knudsen number defined by the length scale 
for variation of the hydrodynamic fields. Note that while 
the strength of the gradients can be controlled by the 
initial or the boundary conditions in the case of elastic 
collisions, the problem is more complicated for granular 
fluids since in some cases (e.g., steady states such as the 
simple shear flow (27l . [28[) there is an intrinsic relation 
between dissipation and some hydrodynamic gradient. 
Here, however we consider situations where the spatial 
gradients are sufficiently small (low Knudsen number). 
Moreover, in ordering the different level of approxima- 
tions in the kinetic equation, one has to characterize the 
magnitude of the external (thermostat) forces relative to 
the gradients as well. As usual, it is assumed that the 
external forces (drag and stochastic forces) do not induce 
any flux in the system and only modify the form of the 
transport coefficients. As a consequence, 7b and ^ are 
taken to be of zeroth order in gradients. 

According to the expansion p6[) for the distribution 
function, the Enskog collision operator and time deriva- 
tive are also expanded in powers of e: 



IV. SMALL SPATIAL PERTURBATIONS 
AROUND THE STEADY HOMOGENEOUS 
STATE 

The steady homogeneous state described in the previ- 
ous Section can be perturbed by small spatial gradients. 
The response of the system to these perturbations gives 
rise to nonzero contributions to the heat and momen- 
tum fluxes, which are characterized by transport coeffi- 
cients. The main goal of this paper is to determine the 
transport coefficients of the driven granular fluid. In or- 
der to obtain them, we consider states that deviate from 
steady homogeneous states by small spatial gradients. 
In this case, the Enskog kinetic equation (0} is solved by 
means of the Chapman-Enskog method |14| conveniently 



JE = jf+ £ 4 1) + ---, d t ^d[ a) +ed\ 1) + ---. (37) 



The coefficients in the time derivative expansion are iden- 
tified by a representation of the fluxes and the cooling 
rate in the macroscopic balance equations as a similar 
series through their definitions as functionals of /. This 
is the usual CE method [H E3] for solving kinetic equa- 
tions. The expansions (|37l) lead to similar expansions 
for the heat and momentum fluxes when substituted into 

Eqs. CEDMSk 



q = q( 0) +eq (1) 



(38) 



In this paper, we shall restrict our calculations to the 
first order in the uniformity parameter e. 
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A. Zeroth-order approximation 

To zeroth order in the expansion, the distribution / 
obeys the kinetic equation 

d (o) .(o) _ 7b u . A m _ 7b _d_ y .(o) 
* to <9V to <9V 

-^^/ (0) = jf [/ (0) ,/ (0) ], (39) 

where J^if^, / (0) ] is given by Eq. (JTSJ) with the re- 
placements n s -> n(r,t) and / s -> /(°)(r,v,t). The con- 
servation laws at this order give 



9 t (0) n = 0. a t (0) u = -^ U; 

to 



n(0) T 



2 ^ - 

7b + TO£ b 

TO 



C (0) T, 



(40) 



(41) 



where is determined by Eq. (|14[1 to zeroth order 
(namely, it is given by Eq. (|31| in the first Sonine ap- 
proximation). The time derivative can be more 
explicitly written as 



dn 

TO 



dUi 1 1 + dT 0t 1 



dV 



T 



(42) 



With this result, Eq. ([39]) becomes 



2 TO o 

-7b - ^Cb + C 
m i 



(0) 



^/ (0) 7b d (0) 
dT to ,9V ' 



^^/ (0) = jf[/ (0) ,/ (0) ]. 



(43) 



It is important to remark that for given values of 7b, 
£b and a, the steady-state condition (f2"U]l establishes a 
mapping between the density and temperature so that 
every density corresponds to one and only one tempera- 
ture. Since the density n(r, t) and temperature T(r, t) are 
specified separately in the local reference state f^°\ the 
collisional cooling is only partially compensated for the 
heat injected in the system by the external driving force 
and so, d^T ^ 0. Consequently, the zeroth-order distri- 
bution function /(°) depends in general on time through 
its dependence on the temperature. On the other hand, 
for simplicity, one could impose the steady-state condi- 
tion (|20|) at any point of the system and so, d^'T = 0. 
This was the choice used in previous theoretical works 
[l5l HH in the case of the stochastic thermostat (7b = 0) 
where the relation = Q°^T was assumed to apply 
also in the inhomogeneous state. As we will see below, 
while the expressions of the shear and bulk viscosities are 



the same in both choices (d ( t 0) T = and d^'T ^ 0), the 
transport coefficients of the heat flux are different. The 
former choice of thermostat (9 t T 7^ 0) will be referred 



3(0) r 



here to as the choice A while the latter (d^T — 0) will 
be referred as to the choice B. Note that the choice A 
has the advantage of a simpler implementation in com- 
puter simulations. However, at the level of kinetic the- 
ory, the fact that d^T ^ gives rise to conceptual and 
practical difficulties not present in the previous analysis 
[T5L [l6j carried out by using the choice B. The above 
difficulties are also present in a recent Chapman-Enskog- 
like method proposed to analyze rheological properties 
around the steady shear flow state [25|, . 

Although the drag parameter 7b and the white noise 
parameter ^ are in general independent parameters, 
henceforth for the sake of simplicity we will assume that 
both parameters are related as 



7b = 13 



(44) 



where {5 is an arbitrary constant. Thus, when /? = our 
thermostat reduces to the usual stochastic thermostat 
[HI while the choice ft = | was previously considered in 
Ref. Q. Here, we will take as arbitrary and only use 
an specific value of j3 at the end of the calculations. 

In the unsteady state, the zeroth-order distribution 
function /W) obeys Eq. (|T9|) . Dimensional analysis re- 
quires that /(°) is also given by the scaled form (|2"2"]l (once 
one uses the relation (1441). namely 



/(°)(r ) v ) t)=n(r ) t)« (r ) t)-V(c,r) 



(45) 



where now c = V /vq, V = v — U being the peculiar 
velocity. Here, the thermal velocity Vo and the reduced 
model parameter £* are defined as in Sec. IIHI with the 
replacement T s — > T(r,t). As in the steady state, the 
temperature dependence of is not only through ^o 
and c but also through £* (see Eq. I|23[)). Thus, 



T 



dT 



il. V/ (0)-V^ (46) 
2 9V J 2 ? d? ' [ ' 



and in dimensionless form Eq. (|39|) can be written as 



1 d 2 

4^ d^^ = Je[<P,<P], 



where 



T 

= -r, C = 



T d-1 



(47) 



(48) 



Upon writing Eq. (|47[) use has been made of the relation 
jSJ. Note that the reduced temperature T* oc £*~ 2 / 3 . 

As before, the explicit form of <p is not known. An 
indirect information on the scaled distribution y> is given 
through its fourth-cumulant 02 (C*) which is defined by 
Eq. (|27p . This cumulant can be obtained by multiplying 
both sides of Eq. (|4"T1) by c 4 and integrating over velocity. 
The result is 
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3d{d + 2) 



[(2(3T* - l)r + Cc 



d(d + 2) 



[Q(l + a 2 )-Ca 2 ] =M4, 



(49) 



r 



where 



I 



dc c e Je[<p, tp\. 



(50) 



In the steady-state, Eq. (j2TJ)) applies and the first term on 
the left hand side of Eq. (|49l) vanishes. In this case, the 
solution to Eq. (j4*§|) is given by Eq. (|5P|) . In general, Eq. 
(|49)l must be solved numerically to get the dependence 
of ai on £* (or equivalently, on the reduced temperature 
T*). An analytical expression of 80,2/8!;* in the steady 
state has been obtained in the Appendix^] Thus, in what 
follows a2(£*) will be considered as a known function of 



V. TRANSPORT COEFFICIENTS 

The analysis to first order in spatial gradients is similar 
to the one worked out in the undriven case 
Some technical details on the determination of the trans- 
port coefficients and the cooling rate are provided in the 
Appendices 151 and ICl The form of the first-order velocity 
distribution function /W is given by 



A(V) • VlnT + B(V) • Vlnn 

+Gi (V) x - Uuj + djUi - ^ y V ■ U 

+V (V)V-U, 



(51) 



where the quantities A(V), S(V), dj (V) and V (V ) 
are the solutions of the linear integral equations (|B13p ~ 
(|B16|) . respectively. However, the evaluation of the trans- 
port coefficients from the above integral equations re- 
quires to know the complete time dependence of the first 
order contributions to the pressure tensor and the heat 
flux vector. This is quite an intricate problem. On 
the other hand, some simplifications occur if attention 
is restricted to linear deviations from the steady state 
described in Section [TTJ In particular, since the kinetic 
and collisional contributions to the heat and momentum 
fluxes are already of first order in the deviations from the 
steady state, one only needs to know the transport coef- 
ficients to zeroth order in the deviations. This means we 
can evaluate the transport coefficients in the steady state 
conditions, namely, when the condition (|20[) applies. 
In this case, Eqs. (|B13|) (IB16|) become 



H 1 



T 

7b 8 
m <9V 



VA 



Z8Q 
2<9£* 
1 
2 



1 



(0) 



7b 



U 



8A 
8V 



F 2 —— 



A+ CA = A, 



(52) 



7b u _ 7b 8 

m <9V m 8~V 

CB = B + C (0) I I 



VB 



8 2 



W 2 
A, 



B 



(53) 



7b u ddj 
m 8Y 



mdV 13 



1 2 8 2 

2 dV 2 = 



7b 



8V 7b 8 



1 



d 2 



8V m 8V 



2 st W 2 



(54) 

V + CV = D, (55) 



where it is understood that the quantities A, B, Cy, and 
D (defined by Eqs. (|B6|) -(|B9 |) . respectively) are evalu- 
ated in the steady state. Consequently, all the transport 
coefficients are given in terms of the steady granular tem- 
perature T s . 

The forms of the collisional contributions to the mo- 
mentum and heat fluxes are exactly the same as those 
previously obtained in the undriven case [I2T [l3| except 
that a 2jS depends on £* (see Eq. (j3"0")) '). Thus, we will 
focus here our attention in the evaluation of the kinetic 
contributions to the transport coefficients and the cool- 
ing rate. Technical details on this calculation are given 
in the Appendix [Cj 



A. Pressure tensor 

To first order in the spatial gradients, the pressure ten- 
sor is given by 

= -V UiUj + 8 3 U t - ^<%V • U) -AJyV-U, (56) 

where 77 is the shear viscosity and A is the bulk viscos- 
ity. While the shear viscosity has kinetic and collisional 
contributions, the bulk viscosity has only a collisional 
contribution. The bulk viscosity A is given by 



where 



2 r(|) 



A-d 



r (d-l)/2 



■\/mT s 



(57) 



(58) 



is the low density value of the shear viscosity in the elastic 
limit. The shear viscosity rj can be written as 



V 



Wo 



2j3m f 2 
2 d-l 



->d~2 



1 



d + 2 



d + 2 



(1 + a)<f>x 



(l + a)(l-3a)0x 



A, 



d + 2 



(59) 
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where h>$ = n s T s /r)o and the collision frequency i/ v is |30| and the value of the derivative (da 2 /d£ t *) s in the steady- 
state is provided in the Appendix [A] Moreover, the col- 



3^o 
Ad 



x [l-a+-d)(l + a)[l + — a 2tS ). (60) 



lision frequency v K is |3C 



1 + a 



V K = Vq- 



d — 1 3 , , „. , , 
— + lg(d + 8)(l-a) 



B. Heat Flux 



296 + 217d- 3(160 + lld)a 



The constitutive form for the heat flux in the Navier- 
Stokes approximation is 



256 



-012,! 



The coefficient /i is 



-kVT — /iVn, 



(61) 



1 + 3^^(1 + a) 



where k is the thermal conductivity and fi is a new coef- where its kinetic contribution is 
ficient not present in the elastic case (a = 1). 
The thermal conductivity k is given by 



(66) 



(67) 



1 + 3^^1 + a) 



AO) _ "^b 



2 2rf+1 (rf- 1) 

(d + 2) 2 7T 



\(i - n) ( l + ^02, s ) /.-o. K»2! 



^Ci 0) (l + Wn X ) + ^-±a 2 . s 
K fo a 

2 d ~ 2 (d - \) / i 

-3 „.V V x(l + «) l + -^ln X 



where 



K 



d(d + 2) 77o 
2(d- 1) m 



d(d + 2) 
a(a — 1)4 



«2,i 



(10 + 2d - 3a + 3a 2 ) j 



(63) 



(68) 



is the thermal conductivity coefficient of an elastic dilute 
gas. The expression of the kinetic part Kk appearing in 
Eq. is 



d-l 



1Kb 

2 T s 



1+30 



M 



/ da 2 



K 



(0) 



C. Cooling rate 

The cooling rate £ is given by 

C = Ci 0) + Cc/V • u. 



(69) 



, 3 / da 2 



. 1 10 + 2d-3a + 3a 2 f da 2 
+a - (1 + a) -4 1 + a C {W 



In Eq. ([Fill), Cs (0) is given by Eq. (|3ljl. 

Cm 



2 d_3 2 At first order in spatial gradients, the proportionality 

^ d + 2 _ ^ constant 6/ is a new transport coefficient for granular 

fluids [HI, [l3j]. For a driven gas, Cu can be written as 



(64) 



where 



16d r(i] 1 a jx ' 



510 



Cf7 — ClO + Cllj 



2 rf - 2 
-3- — <-' 1 



-X0(l-« ), 



(70) 



(71) 



(65) 



9(d + 2)2 d -\ 2 2 . 
Cn = ^ (1 ~ a ) 



32 



2(d+2) 



(1 + a) (I -a) (2a 2 , s - |4* (H?)J 



2m^ o >(0) 



2a 



and the collision frequencies lo and are 

a = (1 + a)i/ {(1 - a 2 )(5a - 1) - ^ [l5a 3 - 3a 2 + 3(4d + 15)a - (20d + 1)] | 

1 + a 



192 



■X^o [30a 3 - 30a 2 + (105 + 24(i)a - 56d - 73] 



(72) 

(73) 
(74) 



I 

Note that the first-order contribution to the cooling trary solid volume fraction <f) and of dilute inelastic gases 
rate vanishes in the limits of elastic gases (a = 1, arbi- 
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(<f> = 0, arbitrary values of the coefficient of restitution 
a) [see Eqs. {nj and ([72])]. 

The expressions for the Navier-Stokes transport coef- 
ficients obtained by using the choice B (i.e., when the 

condition (|TT1) holds locally and so, d^T — 0) are dis- 
played in the Appendix [D] While the expressions of rj 
and A are also given by Eqs. (f57"l) -(f5"9" ]) . the forms of k 
and fx are different to those derived from the choice A. 

In the case = 0, the results derived in this Section 
agree with those obtained in the undriven case [l2l [l3j . 
Another interesting limit corresponds to the low-density 
gas (0 = 0). In this case, A = (jj = while r/, n and fi 
are given, respectively, by 



2/3m C 2 ' 



(75) 



d- 1 



1 + 2a 2 „ 



2 T s 

3 » / 9a2 
2 ?s ^„ 



1 + 3C: 



2C 



(0) 



n s 

( Kk 



3 f r (o) _ !<b 
2 T s 



ci 0) 



f 



-a 2 , 



(77) 



where and ^ K are defined by Eqs. (|60[) and (I66|) . respec- 
tively, with x = I- The expressions (|75|) and (|76|) agree 
with recent results [3ll | derived from the linearized Boltz- 
mann equation for a granular gas heated by the stochastic 
thermostat ((3 = 0). 



VI. COMPARISON WITH COMPUTER 
SIMULATIONS 

The expressions derived in Sec. |V] for the transport 
coefficients and the cooling rate depend on the (steady) 
granular temperature T s , the coefficient of restitution a, 
the solid volume fraction <\> along with the parameter ^ 
characterizing the external energy source. In this Sec- 
tion we will compare our theoretical predictions for the 
thermostats A and B with recent MD simulations carried 
out by Gradenigo et al. 71] for hard disks (d = 2). In 
these simulations, the fluid is also driven by a stochastic 
bath with friction and the two external parameters 7b 
and ^ are related by Eq. (PHI) with P — \- In the steady 
state, they measured the static and dynamic structure 
factors for shear and longitudinal modes for several val- 
ues of the coefficient of restitution a and volume fraction 
<f>. The corresponding best fit of the simulation results 
of the above structure factors allow them to identify the 
kinematic viscosity v = rj/ p, the longitudinal viscosity 



vi = -\ 2— — n + A 

p V d 



(78) 



0.008 



0.006 



0.004 



0.002 




FIG. 3. Plot of the kinematic viscosity v = rj/p as a function 
of the volume fraction cj> for a = 0.6. The solid line is the 
theoretical prediction given by Eq. (|59|) while the dashed line 
is the theoretical result obtained by assuming the elastic form 
of the shear viscosity rj. Symbols are the simulation results 
obtained by Gradenigo et al. 0] from the static (circles) and 
dynamic (triangles) structure factors. 



(76) 



and the thermal diffusivity 



D T = —K. 

an 



(79) 



Figure [3] shows the kinematic viscosity v for disks as a 
function of the volume fraction <j> for a = 0.6. Symbols 
refer to the values of v obtained from MD [7[ by using 
two different procedures (static and dynamical structure 
factors) . As in Fig. [TJ the parameters of the simulation 
are tj, = 1, T(, = 1, m = 1 and a = 0.01. Regarding 
the theoretical results, note that in this case the results 
obtained by using both kind thermostats are the same. 
The theoretical prediction for r\ in the elastic limit (i.e., 
Eq. (J!)!?]) with a = 1 and j h — £ b = 0) but considering the 
a-dependence of the granular temperature given by Eq. 
((32)) is also plotted. This was the theoretical expression 
for v used in Ref. 7] to compare with MD data. At a 
qualitative level, we observe that both theories (the elas- 
tic Enskog theory and the one derived here) reproduce 
in general the trends of simulation data. However, at 
a more quantitative level, it is quite apparent that the 
analytical results obtained in this paper for granular flu- 
ids agree much better with MD than those obtained in 
the elastic case, since the latter clearly overestimates the 
value of v. This is the expected result since the simula- 
tions were carried out for inelastic gases in the presence 
of a stochastic bath with friction. 

The longitudinal viscosity vi is plotted in Fig. Inversus 
the volume fraction (ft for the same systems as in Fig. [3] 
We observe that, in general, the influence of the thermo- 
stat on the longitudinal viscosity is less significant than 
for the kinematic viscosity v since both theories agree 
relatively well. However, the discrepancies with MD are 
more important than in the case of v, specially in the 
low-density limit (<j> = 0.1). While the clastic theory is 
closer to the simulation data than the inelastic theory 
when a = 0.8 (panel (a) of Fig.QJ, the opposite happens 
at a = 0.6 for denser systems (see the panel (b) of Fig. 
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FIG. 4. Plot of the longitudinal viscosity ui as a function of the volume fraction (j> for two values of the coefficient of restitution: 
a — 0.8 (panel a), and a — 0.6 (panel b). The solid lines are the theoretical predictions for i/j obtained by using Eqs. (|57|l and 
(|59|l while the dashed lines are the theoretical results obtained by assuming the elastic forms of the shear viscosity r\ and the 
bulk viscosity A. Symbols are the simulation results obtained by Gradenigo et al. !7] by fitting their numerical data for the 
dynamical correlations of the longitudinal modes. 
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FIG. 5. Plot of the thermal diffusivity Dt = 2n/dn as a function of the volume fraction cj> for two values of the coefficient of 
restitution: a = 0.8 (panel a), and a — 0.6 (panel b). Symbols are the simulation results obtained by Gradenigo et al. [3] 
by fitting their numerical data for the dynamical correlations of the longitudinal modes. The solid lines are the theoretical 
predictions for Dt obtained by using Eqs. (|62[1 - (|64[) . the dotted lines are the theoretical predictions for Dt obtained by using 
Eq. (IDlfl and the dashed lines are the theoretical results obtained by assuming the elastic form of the thermal conductivity k. 




0.5 0.6 0.7 0.8 0.9 1.0 



a 

FIG. 6. (color online) Plot of the of the dimensionless quantity 
71/z/Tk versus the coefficient of restitution a for hard disks 
(d — 2) with a — 0.01 and two different values of the solid 
volume fraction <fr. (a) <f> — 0.1, and (b) <j> = 0.3. The dashed 
line corresponds to the results obtained by considering the 
choice B for <f> = 0.1. Note that jj, = in the elastic case 
(o=l). 



0}. Since the dependence of the shear viscosity 77 on <j) is 
well captured by the inelastic Enskog theory (see Fig. [3]), 
it is evident that the discrepancies between theory and 
MD are essentially due to the bulk viscosity A, whose 
value is specially underestimated at low-density. This is 
a quite surprising result since one would expect that the 
influence of A on the value of v\ increases with increasing 
density since A = for a dilute gas (0 = 0). 

The thermal diffusivity is shown in Fig. [5] for the same 
cases as those considered in Figs. [3] and [4] Surprisingly, 
for strong dissipation and quite dense systems (see the 
panel (b) of Fig. [5j) , the comparison between theory and 
simulation agrees in general better when one uses the 
elastic form for £>t instead of its inelastic expression 
(f55| . These results contrast with the ones recently ob- 
tained 16] for the stochastic driving (i.e., when 7b — > 0, 
keeping 7bTb finite) where it was shown the accuracy of 
the inelastic Enskog theory (see Fig. 1 of Ref. [l6j]) for 
moderate densities and finite collisional dissipation. It 
is important to note that the identification of the trans- 
port coefficients from MD requires to fit the simulation 
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FIG. 7. (color online) Plot of the magnitude of the first- 
order contribution (u to the cooling rate versus the coefficient 
of restitution a for for hard disks (d = 2) with a = 0.01 
and three different values of the solid volume fraction </>'■ ( a ) 
4> = 0.1, (b) <f> = 0.3, and (c) <j> = 0.5. Note that (u = in 
the elastic case (a = 1). 



results for small but not zero values of the wave number 
k. Given that the expressions for the Enskog transport 
coefficients are independent of the wave number (since 
the hydrodynamic regime only strictly holds in the limit 
A; — >■ 0), it is possible that the transport coefficients mea- 
sured in the simulations are still functions of k, specially 
when the smallest value of k considered to get the fit re- 
sults is not close to 0. In particular, the simulation data 
for (j> — 0.3 and 0.5 in the panel (b) of Fig. @] were ob- 
tained for ka = 0.4 and 0.5, respectively. In this sense, 
if one would extrapolate the data shown in Table 3 of 
Ref. 7|, one could conclude that the true value of Dt 
would be smaller than the one shown in this figure when 
ka = 0. More simulations would be needed to clarify this 
point. 

Now we consider the a-dependence of the transport 
coefficient u and the first-order contribution (u to the 
cooling rate. Given that both coefficients vanish in the 
elastic limit, they were also neglected in previous studies 
for heated granular fluids [H, 0] • To assess the impact 
of the term — /zVn in the heat flux, the reduced coeffi- 
cient im/(Tk) is plotted in Fig. [5] versus the coefficient 
of restitution for two different values of the volume frac- 
tion (f> in the case of the choice A. The results derived 
for [i by using the choice B are also plotted for compari- 
son in the case cf> = 0.1. We observe that the coefficient 
/i is negative in the case of the choice B, although its 
magnitude is practically zero. This drawback (/i < 0) of 
choice B is not present in the case of the choice A since 
fj, is always positive for any value of a and </), similarly to 
what happens in the undriven case [HI, EH . In addition, 
although the magnitude of fi is in general smaller than 
that of the thermal conductivity k, we observe that the 
influence of \i on the heat transport could not be con- 
sidered negligible as the degree of dissipation increases. 
The a-dependence of the magnitude of Qu derived from 
the choice A is plotted in Fig. [7] for several values of the 
volume fraction. It is quite apparent that the influence 



of dissipation on | (jj | is more significant than in the case 
of fi, specially at large densities. Consequently, the con- 
tribution of (u to the cooling rate should be considered 
as the rate of dissipation increases. 



VII. LINEAR STABILITY ANALYSIS OF THE 
HYDRODYNAMIC EQUATIONS 

The closed hydrodynamic equations for n, U, and T 
can be obtained by replacing the constitutive forms of 
the pressure tensor, the heat flux, and the cooling rate 
into the balance equations ([T])-©- They are given by 



D t n + nV • U = 0, 



(80) 



7b, 



D t Ui + — U + p^ViP - p _1 V, [77 (ViUj + VjUi 
m 



-6 «V-U 



(81) 



an 



T 



J^V ■ (kVT + pVn) + (Vi£/j + VjUi 



-- % v-u 



X5 i:j V ■ u 



ViUj - TCc/V • U. (82) 



Note that consistency would require to consider up to 
second order in the gradients in the expression (|69[) for 
the cooling rate, since this is the order of the terms in 
Eqs. (|56|) and (|6"T|) coming from the pressure tensor and 
the heat flux, respectively. However, it has been shown 
for a dilute gas that the contributions from the cooling 
rate of second order are negligible [32| as compared with 
the corresponding contributions from Eqs. (|5d| and ([BTI) . 
It is assumed here that the same holds in the dense case 

The form of the Navier-Stokes equations (|50")) - (f52l for 
a driven granular fluid is analogous to that of an ordinary 
fluid, except for the presence of the external bath param- 
eters 7b and ^, the contributions to the cooling rate 
and Cu and the new transport coefficient fi in the energy 
balance equation. In addition, as shown in Sec. [V] and 
depending on the values of the coefficient of restitution 
a, the transport coefficients are in general different from 
those obtained for elastic collisions. 

Equations (|80")) ~ (f52"l) can be linearized around the sta- 
tionary homogeneous state, where the hydrodynamic 
fields take the steady values n s = const., T s = const, 
and U s = 0. A linear stability analysis of the hydrody- 
namic equations (|80p - ([52|) have also been carried out in 
Ref. [7| but neglecting any dependence of the transport 
coefficients on inelasticity and assuming that fi = Cu = 0. 
As mentioned in the Introduction, the only impact of in- 
elasticity on their hydrodynamic equations Q is through 
the cv-dependence of the (steady) granular temperature 
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T s (see Eq. with a 2 , s = 0). Thus, it is worth to 
assess to what extent the previous theoretical results [7[ 
are indicative of what happens when the correct expres- 
sions for the transport coefficients and the cooling rate 
are considered. This is the main motivation of this Sec- 
tion. 

We assume that the deviations Sy a (r,t) = y a (r,t) — 
y sa {t) are small, where Sy a (r, t) denotes the deviations of 
n, U, and T from their values in the steady homogeneous 
state. To recover previous linear stability results [33[ 
derived in the undriven case, let us consider the following 
(reduced) time and space variables: 



1 

T = -U s n 



(83) 



The dimensionless time scale r is a measure of the aver- 
age number of collisions per particle in the time interval 
between and t. The unit length introduced in the sec- 
ond equality of (1531 corresponds to the mean free path 
of gas particles. 

A set of Fourier transformed dimensionless variables 
are then introduced by 

( \ 6n ^ ( \ M a ( ^ gW 
Pk(r) = ~n ' Wk(r) = Fr~l — ' 0k(r) = ~^T~> 

(84) 

where 5y ka = {p^, Wk(r), 0u.(t)} is defined as 

■ r '5y a (r',r). (85) 



Sykair) = J dr' 



Note that in Eq. (|85|) the wave vector k is dimensionless. 

In Fourier space, as expected, Eq. (|5Tj) shows that the 
d — 1 transverse velocity components Wkj_ = Wk — (wk • 
k)k (orthogonal to the wave vector k) decouple from the 



other three modes and hence can be obtained more easily. 
Their evolution equation can be written as 



A + 2 7 * + I»,**c 2 ) w k± 0. (Mi) 



where 

* Tb * V {0.7^ 

The solution to Eq. ([551) is 

w k ±(k, r) = w k± (0) cxp [A±(k)r] , 

where 

1 



A j/,; : (•-'- • • -r,*k 2 



(89) 



Since 7* and 77* are both positive, then A±(k) becomes 
negative for any finite wave number k and so the transver- 
sal shear modes of the driven gas are linearly stable. This 
result contrasts with the ones obtained in the undriven 
case [33[ where it was shown that the transversal shear 
modes become unstable for values of k smaller than a 
certain critical wave number. 

The remaining (longitudinal) modes correspond to pk, 
6*k, and the longitudinal velocity component of the veloc- 
ity field, u>k|| = Wk • k (parallel to k). These modes are 
coupled and obey the equation 



dSy ka (T) 

— = M a0 Sy kp {T), 



(90) 



where <5j/ kQ (r) denotes now the set {p^, 6^, w^n } and 
is the square matrix 



M 



M 





—ikp*C p 



-3Co 





■4 7 * - 
—ikp* 



-ik 



D* T k 2 



-2 7 * - 



(91) 



where Here, p s = mn s is the mass density. In the above equa- 

(0) tions, it is understood that the transport coefficients 77, 

Q = ^ s | (92) vu Dt : and p are evaluated in the homogeneous steady 

^/2T s /m state. In addition, the quantities g{4>) and C p (a,4>) ap- 

pearing in the matrix M are given by 



n s l s 

and 



l + 2 a - 2 {l + a) X <t>, (93) q 

9(0) = l + ^lnxW, (96) 



PsVl r>* - nsL>T (qa\ 



2a 1 d yfthT s 2a 1 - d y /T s /m Cp (a, (j,) = 1 + 6^7 lnp*(a, <\>) 

dtp 

"' =M- (95) = i + ^)____^L_^ ) (97) 



da l - d T s ^/^T s ^' K ' \ + 2 d - 2 {\ + a)4>xW 
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where in the last equality use has been made of the ex- 
plicit expression of p* given by Eq. (|93p . If one assumes 
fx* = = 0, the matrix (|91|) agrees with the dynamical 
matrix obtained when the gas is heated by a stochastic 
thermostat (7b = but jbTb = finite) [5(. 

The longitudinal three modes have the form 
exp[Af(fc)r] for £ — 1,2,3, where A^(fc) are the eigen- 
values of the matrix M , namely, they are the solutions of 
the cubic equation 



where 



A 3 + A(k)A 2 + B(k)A + C(k) = 0, 



A(k) = 3(C * + 2 7 *) + fc 2 (!/? + £>?.), 



B(k) = k 4 ^D* T + k 2 

+2y*D^ + (3Co + 4 7 *)^ 



P*C P + P * [ ~ d P*+Cu 



(98) 



(99) 



(100) 



C{k) = p*k 2 [C p (4 7 * + 3Co) - 2 ff Co + (C P D* T - //) k 2 ] . 

(101) 

One of the longitudinal modes (the heat mode) could be 
unstable for k < fch, where fch is obtained from Eq. (|98|) 
when s = 0, namely, C(kh) = 0. The result is 



k 2 



(102) 



On the other hand, an analysis of the dependence of fc 2 , on 
the coefficient of restitution a and the volume fraction 4> 
shows that fc 2 , < for any value of a and <f>. Thus, there 
are no physical values of fch for which the heat mode 
becomes unstable. Consequently, all the eigenvalues of 
the dynamical matrix M have a positive real part and no 
instabilities are found due to the presence of the external 
bath. This conclusion agrees with the results obtained in 
Refs. [3 and for driven granular fluids. 

In summary, the results obtained here including the 
complete a-dependence of the transport coefficients show 
no new surprises relative to the earlier works pi @] by 
considering the elastic Enskog expressions for the above 
coefficients. Of course, the quantitative forms for the dis- 
persion relations can be quite different in both (elastic 
and inelastic) approaches since the impact of dissipation 
on the transport coefficients and the cooling rate is sig- 
nificant and so, their functional forms differ appreciably 
from their elastic forms. 



VIII. DISCUSSION 

In this paper, we have determined the transport co- 
efficients of a granular fluid driven by a stochastic bath 
with friction. The results have been obtained within the 
framework of the (inelastic) Enskog kinetic theory and 
they are expected to apply over a wide range of densi- 
ties. Our goal is not only academic since, from a prac- 
tical standpoint, many of the simulations reported [3-@] 



for flowing granular materials have used external driving 
forces to fluidize the system. For this reason, it would 
be convenient to provide to simulators with the corre- 
sponding expressions of the transport coefficients when 
the granular fluid is heated by a thermostat. In fact, due 
to the lack of the above expressions, in most of the cases 
it is assumed that the forms of the transport coefficients 
of the driven granular fluid are the same as those given 
by the elastic Enskog theory (jyj. However, as expected 
from previous theoretical works [TH, HU , the present re- 
sults show again that the expressions for the transport 
coefficients clearly differ from those obtained for ordinary 
fluids so that, one should use the true inelastic Enskog co- 
efficients to analyze granular flows driven by thermostats. 

The transport processes considered are those for a 
driven fluid with small spatial gradients of the hydrody- 
namic fields. In this situation, the Enskog equation has 
been solved by means of the Chapman-Enskog method 
[lij up to the first order in the spatial gradients. Since 
these gradients have been assumed to be independent 
of the coefficient of restitution, although the correspond- 
ing hydrodynamic equations restrict their applicability to 
first order in gradients, the transport coefficients appear- 
ing in these equations hold a priori to arbitrary degree 
of dissipation. 

An important but subtle point is the generalization of 
the driving external forces (which are mainly used in ho- 
mogeneous situations) to inhomogeneous states. This is a 
crucial step since one has to consider situations close to 
steady homogeneous states to determine the transport 
coefficients from the Chapman-Enskog expansion. Al- 
though the above generalization is a matter of choice, it 
has important implications in the final expressions of the 
transport coefficients. For simplicity, in previous works 
on heated granular gases [HI, [l6[ it was assumed that 
the external driving force has the same form as in the 
homogeneous case, except that their parameters are lo- 
cal quantities. As a consequence, the parameters of the 
force are chosen to impose a stationary temperature in 
the zeroth-order solution (i.e., df'T = 0). However, for 
general small perturbations around the steady homoge- 
neous state, it is expected that the density and temper- 
ature are specified separately in the local reference state 
f(°' and so, the temperature cannot be stationary at any 
point of the system (i.e., d^T ^ 0). This choice is more 

general than the previous one (d^T — 0) and has the 
advantage of a simpler implementation on computer sim- 
ulations since the parameters of the driven external force 
are constant, even for inhomogeneous states. 

As mentioned in the Introduction, the fact that 
d[ 0) T £ gives rise to conceptual and practical difficul- 
ties not present in the case of the choice B. One of them 
is that the evaluation of the complete nonlinear depen- 
dence of the transport coefficients on dissipation requires 
in principle the analysis of the hydrodynamic behavior 
of the unsteady reference state. This involves the corre- 
sponding numerical integration of the differential equa- 
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tions obeying the velocity moments of the zeroth-order 
distribution /(°) (see for instance, Eq. (l4"9l for the fourth 
degree moment a-i of Z*- '). This is quite an intricate long 
problem. However, given that here we are interested in 
evaluating the momentum and heat fluxes in the first or- 
der of the deviations from the steady reference state, the 
transport coefficients must be determined to zeroth or- 
der in the deviations. As a consequence, the steady-state 
condition (|20|) applies and the transport coefficients and 
the cooling rate can be defined in terms of the hydro- 
dynamic fields in the steady state. Explicit expressions 
for these quantities have been obtained after consider- 
ing the leading terms in a Sonine polynomial expansion. 
These explicit forms have been displayed in Sec. [V] and 
Appendix ID1 for the choices A and B, respectively. More 
specifically, in the case of the choice A, the bulk A and 
shear 77 viscosities are given by Eqs. (f57|) and (|59|) . respec- 
tively, the thermal conductivity k is given by Eqs. (|62[) 
and the coefficient \i is given by Eqs. (|57| and 
and the cooling rate £ is defined by Eqs. (j6"9"j) - (j?4")) . All 
these expressions clearly show the complex dependence of 
the set {A, rj, k, (a, £} on the granular temperature T, the 
coefficient of restitution a, the solid volume fraction <p 
and the model parameter In the case of the choice B, 
our results show that the expressions of A and 77 are the 
same as those obtained from the choice A but the forms 
of k and fi are different (they are given by Eqs. (|D1[) and 
(|D2[) , respectively). An important drawback of the re- 
sults derived from the choice B is that the coefficient /_t 
can be negative (see Fig. [SJ) , although its magnitude is 
very small. 

A comparison with recent MD simulations Q carried 
out for a granular fluid driven also by a stochastic bath 
with friction has been made in Sec. IVI1 The comparison 
has been displayed in Fig.[3]for the kinematic viscosity v, 
Fig. |U for the longitudinal viscosity vi and Fig. [5] for the 
thermal diffusivity Dt- It is quite apparent that while 
the predictions of the driven kinetic theory compares very 
well with simulation data for v in a wide range of den- 
sities, some discrepancies appear in the cases of v\ and 
Dt as the gas becomes denser. Surprisingly, in the case 
of Dt, the comparison agrees better when one uses the 
elastic form of Dt in the more inelastic system (a = 0.6) 
studied. We think that this disagreement is in part due 
to the fact that while the simulation data have been ob- 
tained for small but finite values of the wave number k, 
the Enskog expressions for the transport coefficients only 
strictly apply in the limit k — > 0. Moreover, given that 
these discrepancies appear at sufficiently high densities, 
it could also reflect the limitations of the Enskog equa- 
tion (which is based on the molecular chaos hypothesis) 
as the granular fluid becomes denser. 

With these new expressions for the momentum and 
heat fluxes and the cooling rate, a closed set of hydro- 
dynamic equations for situations close to homogeneous 
steady states has been derived. A stability analysis of 
these linearized hydrodynamic equations with respect to 
the homogeneous steady state has been carried out to 



identify the conditions for stability in terms of dissipa- 
tion. Our results show that the driven homogeneous state 
is stable for any value of dissipation at sufficiently long 
wavelengths. This conclusion agrees with previous find- 
ings [H,[3] obtained by using the elastic expressions of the 
transport coefficients. 

Finally, it would be interesting to extend the results 
derived in this paper for a single granular gas to the im- 
portant subject of granular mixtures. Given the difficul- 
ties associated with multicomponent systems, the tracer 
diffusion could be perhaps a good starting point to pro- 
vide some insight into the general problem. Work along 
this line will be carried out in the near future. 
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Appendix A: Behavior of the fourth-cumulant of the 
zeroth-order distribution in the vicinity of the 
steady state 

Although the determination of 02 (^*) requires numeri- 
cal work, one can obtain analytically this quantity in the 
vicinity of the steady state by means of the derivative 
da2/d^* evaluated at the steady state. This derivative 
appears in the expressions of the thermal conductivity k 
(see Eq. (JM])) and the first-order contribution £[/ to the 
cooling rate (see Eq. (|72"jl ). This Appendix addresses the 
evaluation of the above derivative. 

In order to determine da2/d^* from Eq. (j4"9")l . we first 
assume that <p can be well described by the lowest So- 
nine approximation (|29|) . Then, approximate forms for 
Co = (2/d)n2 and LI4 can be obtained when one uses the 
distribution (129p and neglects nonlinear terms in a^- The 
results are [17j 



,, , „(0) , „ , ,,(0) (1) 



where 



u (0) 



V2r(f 



-x(i-« 2 ), 



/4 1} 



3 (o) 
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(A3) 
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(10d + 39 + 10a 2 ) + 



1 
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(A4) 



With the use of the approximations (IA1[) and retaining 
only linear terms in a 2 in Eq. (|49p . the derivative da 2 /d^* 
is given by 



da 2 



M i 0) -(<i + 2) M W- 


j§(d + 2)/4°> 


(1) d{d+2) c * 
2 s 


02 


3d(d+2) tt 
8 ? 


(2^t* - i)e + 1 


i M 2 + a 2y 







(A5) 



However, some care must be taken in Eq. (|A5[) at the 
steady state since the numerator and denominator of Eq. 
(|A5[) vanish and so, the corresponding expression for the 
derivative da 2 / becomes indeterminate. This diffi- 
culty can be solved by means of l'Hopital's rule. After 
some algebra, it is straightforward to see that the steady- 
state value of the derivative A = (da 2 /d^*) s obeys the 
quadratic equation 



|(d + 2) M «CA 2 



d(d+2) 



A- 



(l + 2/?T s *K, 
d(d + 2 



a 2 , s = 0, (A6) 



where T* = T s /T\,. Since a 2 . s is in general very small, it 
is expected that the magnitude of A be also quite small. 
An analysis of the solutions to Eq. (|A6I) shows that in 
general one of its roots is much larger than a 2jS while the 
other is of the order of a 2 , s - We take the latter one as 
the physical root of the quadratic equation (IA6|) . 



Appendix B: First order approximation 

The application of the Chapman-Enskog method up 
to the first order approximation follows similar mathe- 
matical steps as those made before in the undriven case 
HI • Some details on this derivation are provided 
in this Appendix. Up to the first order in the expansion, 
the velocity distribution function /W obeys the kinetic 
equation 



-f 2 



d 2 



777 <9V 777 9V 



dV 2 



/ (1) = -(a t (1) + v.v) /<°>-j£>[/]. 



(Bl) 

Here, J$p[f] means the first order contribution to the 
expansion of the Enskog collision operator and C is the 
linear operator 



(i) 



(jf [/ (0) ,/ (1) ] + jf [/ (1) ,/ (0) ] 



(B2) 



The macroscopic balance equations to first order in the 
gradients are 



D[ 1] n = -77V ■ U, Df ) U l = -(mn) _1 ViP, (B3) 



.BE 
dn 



V-U-C (1) T. 



(B4) 



where Df' = df' + U ■ V and £W is the first order 
contribution to the cooling rate. Use of Eqs. (|B3[) in (|B1[) 
and taking into account the form of [/] obtained in 
Ref. [l2j for the undriven case, one gets 

^(o) + £ ) / (i ) _2ku.^-2k » . v/w 

V / 777 OV 777 W 

-^b^2/ (1) =A-VlnT + B-Vln77 

+Cii^ (ftt/j + SjUi - • u) +DV ■ U, (B5) 

where 

A(V) = -VTd T f {0) - -V v / (0) - /C [Tdr/ (0) 1 , (B6) 
P L J 



B (V) = -V/W - ^ (l + ^ InP* ) Vv/ (0) ) 



1 



/( o) 
A; 

ay/ 



(0) 



(B7) 
(B8) 



Here, V v = d/dV, 



= l + 2 d - 2 (l + a) X 0, 



TlT 



(B9) 
(BIO) 

4> is defined by Eq. (|33|) and JCi is the operator 

K,i[X] = a d x y" dv 2 f daQ(a ■ gia)(ff ■ gw)^ 

x a- 2 / (0) (vi')X(v^) + / (0) (vi)X(v 2 )' , 

(Bll) 

where v" = vi — |(1 + a -1 )^ • gia)? and v 2 ' = V2 + 
i(l + a -1 ) (3 s ■ gi 2 )3\ In Eq. ([B9]). Cc/ is defined through 
the expression 



C = C (0) + Cc/V • u, 



(B12) 
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where is the cooling rate evaluated at zeroth-order. 

The solution to Eq. (|B5|) can be written in the form 
(|5ip . where A, 13, Cij, and T> are unknown functions 
of the peculiar velocity. Since the gradients of the hy- 



drodynamic fields are all independent, Eq. (|B5[) can be 
separated into independent equations for each coefficient. 
This yields the following set of linear, inhomogeneous in- 
tegral equations: 



ml j to ox 



^fi-~8f) + k (0) 



T 



dV 



\ m T J to ox to ox 2 



d 2 



CB = B + C (0) I 1 + <j>i- lux ) A (B14) 
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27b _ rn^b 
to T 



+ C (0) ) TdrD 
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m dV 



2 st W 2 " 
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Upon deriving Eqs. (|B13|) — (|B16|) use has been made of 
the result 



<9 t (0) VlnT = Vd y t'> hxT = V 



(o) 



\ T m J 



-C (0) (l + ^lnx)Vlnn 
- C 2 fl_3agVl c( o) 



T 



2 8C 



VlnT. 
(B17) 



As noted in Section [V] in the first order of the devia- 
tions from the steady state, we only needs to know the 
transport coefficients to zeroth order in the deviations 
(steady state conditions). This means that the term 

27b _ m£l + (0) 
to T 



appearing in the left-hand side of Eqs. (IB13[) - (IB16[) van- 
ishes. The differential equations for the transport coef- 
ficients thus become simple coupled algebraic equations. 
They are given by Eqs. ([52]) - ([551) . 



Appendix C: Kinetic contributions to the transport 
coefficients 



In this Appendix we determine from Eqs. ()52j) -(|55 j) the 
kinetic contributions to the transport coefficients to n and 
\x as well as the first order contribution (u to the cooling 
rate. Given that all these coefficients are evaluated in the 
steady state, the subscript s appearing along the main 
text will be omitted in this Appendix for the sake of 
brevity. 



r 



We start with the shear viscosity 77. Its kinetic part ijk 
is given by 

^ = - {d -i)( d+ 2) J *^CV), (CD 



where 



D tj = to V,\{, 



(C2) 



To obtain 77^, we multiply Eq. (|54[) by Dij and integrate 
over velocity. The result is 



27b 



1 



(d-i)(d + 2) 
d 



f 



(0) 



j dVAi(v) 

(C3) 



where 



(C4) 



J dvAj(V)£C -(V) 



The collision integral of the right hand side of Eq. (|C3 
has been evaluated in previous works [IH, [l3j : 



dV/ 



= 2 d - 2 (d-l)nT X <P 
x (l + a)(l - 3a). (C5) 
Thus, the kinetic part rjk can be written as 



nT 



27b 



id-2 



d+2 



(l + a)(l-3a)^x 



(C6) 



where use has been made of the relation (PPT| . In order 
to get an explicit expression for 77^, one has to evaluate 
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the (reduced) collision frequency v* . It can be evaluated 
by considering the leading terms in a Sonine polynomial 
expansion of the function (V). Here, we have consid- 
ered a recent modified version of the standard method 
[30l HH that yields good agreement with computer sim- 
ulations even for quite strong values of dissipation [361 ]. 
The final form of the shear viscosity 77 is given by Eq. (|59|) 
when one takes into account the steady-state condition 
(®. 

The kinetic parts Kk and fik of the transport coeffi- 
cients characterizing the heat flux are defined, respec- 
tively, as 

K k = ~ f dvS(V)-A(V), (C7) 



«fe = -^/ dvSCV)-B(V), 



where 



S(V) = ( —V - 



in rr2 d- 



T V. 



(C8) 



(C9) 



We obtain first the kinetic part Kk- It is obtained by 
multiplying Eq. ([521 by S(V) and integrating over V. 
The result is 



lm 3 (, , M 



2 t v 1 + 3 9e 



2C 
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dVS(V) • A, 
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where 



f dvS(V) ■ £A(V) 
JdvS(V)A(V) ■ 



(Gil) 



The right hand side of Eq. (|C10|) is given by 



-I?!**** j dVS(V).K[/(°) 
2^ d£* da 2 J l 



dVS(V) • K 



'v • 



(v/ 
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where use has been made of Eq. (prol) . The last two terms 
on the right hand side of Eq. (IC12[) can be evaluated more 
explicitly and the result is [29( 

nT 2 



where kq is the low density value of the thermal conduc- 
tivity of an elastic gas (defined by Eq. (l63l ) and 
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The expression (|C15|) for Kk is still exact. In order to 
get explicit results, one considers the form (|3"Tj) for £W 
3q,2-j and evaluates v K by considering again the leading terms 
-I in a Sonine polynomial expansion of A(V). With these 
(C14) approaches, one gets the expression (|6"6"|) for while 



With the above results, the kinetic part Kk can be finally 
written as 

-1 



1 



K k = KoVq- 

14 

2a 



+ i^fl + 3^° 
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where Cm is defined by Eq. ([55]). Use of (|C17|) in Eq. 
(|C15|) gives the final result. 



To evaluate multiply Eq. (|53[) by S(V) and inte- 
(C15)grate over velocity to get 
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The last term on the right hand side of Eq. (|C19|) is 
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The final expression of /Zfc is obtained from Eq. (|C19[) 
when one substitutes Eq. (|C14[) into Eq. (|C19|) . How- 
ever, this expression is not explicit unless one knows the 
collision frequency v^. To determine it, one takes the 
leading terms in a Sonine polynomial expansion of B(V) 
and gets v >i = v K . This finally yields Eq. (1551) for jUfc. 

We consider finally the first-order contribution (u to 
the cooling rate. This coefficient is given by Eq. ([70| . 
where 



1 



,(d-l)/2 



Xm(l — a 2 ) 



Cn = 



and the unknown T> verifies the integral equation (|55p. 
An approximate solution to (|55p can be obtained by tak- 
ing the Sonine approximation 
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Substitution of Eq. fC22|) into Eq. (fC2"T|) gives 
3(d + 2) . 2 . / 3 \ 
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The coefficient is determined by substituting Eq. 
(|C22[) into the integral equation (|55)l . multiplying by 



F(V) and integrating over V. After some algebra one 
gets the expression ([72")) for £n. 



Appendix D: Expressions for choice B 

In this Appendix we display the expressions for the 
Navier-Stokes transport coefficients ?y, A, k, and fi by 
using the choice B defined by the condition d^T = 0. 
The application of the Chapman-Enskog method to this 
case follows similar mathematical steps as those made 
for the choice A (<9 t (0) T ^ 0). The results show that the 
expressions of n and A are the same as those obtained for 
the choice A [see Eqs. ([57)1 - ([59")) ]. However, the forms of 
k and u are different since they are given by Eqs. (|62[) 
and ([S7]) . respectively, but their corresponding kinetic 
contributions are 
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where the collision frequency v K is defined by Eq. (|66|) . 
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